Prictionless atom cooling in harmonic traps: 
a time-optimal approach 
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In this article we formulate frictionless atom cooling in harmonic traps as a time-optimal con- 
trol problem, permitting imaginary values of the trap frequency for trasient time intervals during 
which the trap becomes an expulsive parabolic potential. We show that the minimum time so- 
lution has "bang-bang" form, where the frequency jumps suddenly at certain instants and then 
remains constant, and calculate estimates of the minimum cooling time for various numbers of such 
jumps. A numerical optimization method based on pseudospectral approximations is used to obtain 
suboptimal realistic solutions without discontinuities, which may be implemented experimentally. 
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I. INTRODUCTION 

Frictionless atom cooling in a harmonic trap is defined 
as the problem of changing the harmonic frequency of the 
trap to some lower final value, while keeping the popu- 
lations of the initial and final levels invariant, thus with- 
out generating friction and heating. Achieving this goal 
in minimum time has many important potential applica- 
tions. For example, it can be used to reach extremely low 
temperatures inaccessible by standard cooling techniques 
[l[ , to reduce the velocity dispersion and coUisional shifts 
for spectroscopy and atomic clocks and in adiabatic 
quantum computation 13;]. It is also closely related to 
the problem of moving in minimum time a system be- 
tween two thermal states, as for example in the transition 
from graphite to diamond Q. By using optimal control 
theory, it was proved that minimum transfer time can 
be achieved with "bang-bang" real frequency processes, 
where the frequencies change suddenly at certain instants 
and then stay constant In another recent paper 
it was shown that when the restriction for real frequen- 
cies is relaxed, allowing the trap to become an expulsive 
parabolic potential at some time intervals, shorter trans- 
fer times can be obtained. Based on the theory presented 
in 0, in this article we reformulate the frictionless cool- 
ing problem as a minimum-time optimal control problem, 
permitting the frequency to take real and imaginary val- 
ues in specified ranges. We then show that the optimal 
solution has "bang-bang" form and use this fact to calcu- 
late estimates of the minimum transfer time for various 
numbers of frequency jumps. We finally use a numerical 
optimization method based on pseudospectral approxi- 
mations to find suboptimal realistic solutions which do 
not suffer from discontinuities and are thus appropriate 
for experimental implementation. The efficiency of the 
method is demonstrated by several numerical examples. 



II. FORMULATION OF THE COOLING 
PROBLEM IN TERMS OF OPTIMAL CONTROL 

Consider the one-dimensional time-dependent har- 
monic oscillator with Hamiltonian 



1 ^9 •muP'{t) ^2 



(1) 



with initial frequency w(0) — loq a.t t — Q and final fre- 
quency u}(tf) = ujf < itJo at the final time tf. This cor- 
responds to a temperature reduction by a factor ujf /ujq. 
The goal is to find a path Lo{t) between these two val- 
ues so that the populations of all the oscillator levels 
n = Q, 1,2, ... dX t = tf are equal to the ones at t = 0. 
We would also like to achieve this in minimum time tf. 
It was shown in 5] that appropriate w(t) can be effi- 
ciently engineered by using an invariant of the motion 
([1]). Additionally, by relaxing the restriction > 0, 

allowing {t) < for some time intervals where the po- 
tential becomes expulsive, shorter cooling times can be 
obtained. In the following we present an overview of the 
corresponding theory, which will lead naturally to the 
formulation of the problem in terms of optimal control. 
The basis of the analysis is the invariant of the motion 
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where tt — bp — mbq plays the role of a momentum conju- 
gate to q/b and the dots represent derivatives with re- 
spect to time 0. The scaling dimensionless function 
b — b{t) satisfies the subsidiary condition 



b + uj^{t)b. 
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an Ermakov equation where real solutions must be cho- 
sen to make / Hermitian. I{t) has the structure of a 
harmonic oscillator Hamiltonian, with time-dependent 
eigenvectors \n{t)) and time- independent eigenvalues 
(n -I- l/2)hLUQ. The general solution of the Schrodinger 
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equation is a superposition of orthonormal "expanding 
modes" 



(4) 



where a„(i) = — (n + l/2)aJo /J dt' /ti^, and c„ are time- 
independent amplitudes. For a single mode, 
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where is the Hermite polynomial of degree n. The 
time-dependent average energy of the mode is 



(2n + l)h 
4^0 
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The average position is zero and the standard devia- 
tion a = {J x^l'i'nl'^dxy^^ is proportional to 5, cr = 
b{n + l/2y/^/{mujQ/hy/^, which underlines the physi- 
cal meaning of the scaling factor. 

The approach taken in Q is to leave u}{t) undetermined 
at first and impose properties on b and its derivatives at 
the boundaries t = and t = tf to assure that: 

1. Any eigenstate of H{0) evolves as a single expand- 
ing mode 

2. This expanding mode becomes, up to a position- 
independent phase factor, equal to the correspond- 
ing eigenstate of the Hamiltonian H{tf) of the final 
trap. 

When the above are satisfied, the populations in the in- 
stantaneous basis are kept equal at the initial and final 
times. It is not hard to find the corresponding bound- 
ary conditions. By choosing b{0) = 1, 6(0) = at i = 0, 
H{0) and J(0) commute and have common eigenfunc- 
tions at that instant. Since uj{0) = ojo, it holds that 
6(0) — from ([3]). These conditions imply that any initial 
eigenstate of H{0) will evolve according to the expanding 
mode ([5|). In general H{t) and I{t) will not commute for 
t > 0. At t = tf it is desirable for 'i'nitfjx) to be pro- 
portional, up to the global phase factor e*""^*^-*, to the 
corresponding eigenstate of the final trap. If we impose 
b{tf) 7 = (wo/w/)i/^6(t/) = 0,b{tf) = 0, then from 
([3]) we get io{tf) — ujf and from ([5]) we see that 4'n(i/, x) 
has the desired form. After fixing b(t) and its derivatives 
at the boundaries, b{t) can be chosen as a real function 
satisfying these conditions. For example, substituting the 
simple polynomial ansatz 



into the six boundary conditions gives six equations that 
can be solved to provide the coefficients, 

b{t) = 6(7 - 1)5^ - 15(7 - !)«'' + 10(7 - l)s^ + 1, (8) 

where s — t/tf. Once b{t) has been determined, the 
physical frequency Lu(t) is obtained from the subsidiary 
condition ([3]). 

Note that in the above method, the duration i/ is con- 
sidered to be fixed and there are no bounds on the fre- 
quency uj{t). An alternative approach is to express the 
frictionless cooling problem as a minimum-time optimal 
control problem, incorporating possible restrictions on 
uj(t) due for example to experimental limitations. If we 
set 



Xi = 6, X2 



UJo 



-, u{t) 



(9) 



and rescale time according to inew = "^o^old' obtain 
the following system of first order differential equations, 
equivalent to the Ermakov equation (jS]) 



Xi = X2, 
X2 = —UXi 



(10) 
(11) 



The optimal control problem is: Find —vi < u{t) < V2 
with m(0) = l,u{tf) = 1/7^ such that starting from 
(xi(0),X2(0)) = (1,0), the above system reaches the fi- 
nal point (a;i(i/),X2(i/)) = (7,0) in minimum time tf 
(note that 7 = (wo/w/)^/^ > xhe boundary con- 
ditions on the state variables {xi,X2) are equivalent to 
those for b and b, while the boundary conditions on the 
control variable u lead to the corresponding conditions 
for b. Parameters f 1 , f 2 > define the allowable values 
of u{t) with V2 > m(0) = 1. Note that the possibility 
a;2(t) < (expulsive parabolic potential) for some time 
intervals is permitted. Finally observe that the above 
system describes the one-dimensional newtonian motion 
of a unit-mass particle, with position coordinate xi and 
velocity X2- The acceleration (force) acting on the parti- 
cle is —uxi + l/ccf . This point of view can provide useful 
physical insight, as we will see later. 

The advantage of expressing the cooling problem in 
terms of optimal control is that analytical and numerical 
tools from this area can be used to engineer uj{t), while 
taking into account possible limitations on the frequency. 
The control-theoretical framework has been successfully 
employed to solve various problems in quantum dynamics 
7-25]. We show how this can be done for the problem 
at hand in the following sections. 



III. THEORETICAL OPTIMAL SOLUTION OF 
BANG-BANG TYPE 



6(i) = Ea,t^" 



j=0 



(7) The form of the theoretical time-optimal solution can 

be found using Pontryagin's maximum principle |26j . 
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which we state here in order to keep the paper self- 
contained. 

Maximum Principle for Time-Optimal Prob- 
lems: Consider the autonomous dynamical system 



X = f(x, u), 



(12) 



where x = {xi,X2, . ■ ■ ,Xn) G X (state space), u = 
it2, . . . , Um) S U (control region), and f = 
(/i, /2, • ■ • , fn), with functions /i(x, u) continuous in the 
variables x, u and continuously differentiable with re- 
spect to X. The corresponding control Hamiltonian is 
defined as 



Hcip,x, u) = ^p„/„(x, u), 



(13) 



where p = {pi,p2, ■ ■ ■ ,Pn) is the adjoint vector. Let 
u(t),0 < t < tf, he a.n admissible control which trans- 
fers the state vector from xg to x/, and let x{t) be the 
corresponding trajectory, so that x(0) = xo,x(t/) — 
Xf. For u(t),x(t) to be time-optimal, it is necessary 
that there exists a nonzero, continuous vector function 
p(i) = {pi{t),p2{t), . . . ,Pn{t)) such that: 



1. 



dp 
9x 



(14) 
(15) 



The first equation is equivalent to the system equa- 
tion ([12]), while the other is the equation for the 
adjoint vector. 

2. For all < t < the function iJc(p(i), x(t), u) 
of the variable u e U attains its maximum at the 
point u = u(i). 

3. Hc{p{t),'x.{t),\i{t)) = c > 0, c constant. 

For the system described by (fTOj) and ([Tl]) the states 
{xi,X2) € X = (0, +CX)) X R and the control u G U = 
[— i;i,W2]. Note that xi{t) > because the starting point 
is a;i(0) = 1 > 0, the evolution is continuous for xi ^ 
and when a;i — >■ 0+ there is a "repulsive force" ~ l/x^ 
that forces xi to increase. The system equations satisfy 
the necessary smoothness conditions in spaces X, U. The 
control Hamiltonian is 



Hc{pi,P2,Xl,X2,u) = piX2 

Substituting (IT6t into (fT5)) gives 
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(17) 
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According to the maximum principle, point 2 above, 
the time-optimal control u{t) maximizes the control 



Hamiltonian at each time. Note that He is a linear 
function of the control variable u. Since u is bounded, 
—vi<u<V2, the optimal control that maximizes He 
is determined by the sign of the coefficient of w, which 
is —P2X1. But xi > 0, thus when p2 0, the optimal 
control in {0,tf) is given by 



When 



-Vi, P2> 
V2, P2 <0 



P2=0 



(19) 



(20) 



for some time interval, the maximum principle provides 
a priori no information about the optimal control in this 
interval, which in that case is called a singular control. 
In general, singular extremals can play some role in the 
control of quantum systems [l^, . We show that this 
is not the case for our problem, i.e. that condition (PO)) 
cannot hold for any time interval [^1,^2] C (0,^/). Sup- 
pose that P2{t) — ioi t G [^1,^2], then from (fT8|) we 
have pi — ~p2 = for i g [ii, ^2]- Thus pi{t) — P2{t) = 
for t G [^1,^2], in contradiction with the maximum prin- 
ciple that requires the vector p(i) — {pi{t),p2{t)) to be 
nonzero. So p2 can be zero only at specific moments 
(switching times). The optimal control has "bang-bang" 
form (|19p , where the controller changes from one bound- 
ary value to the other at the switching times. 

Observe that when u is a constant and Eqs. ([TOl) and 
PT|) are satisfied, then 



1 



(21) 



where c is a constant. From ([6]), ([9]) and ([2T|) we find that 
{H{t))n/fvjJo = (2n -t- l)c/4, so the paths of constant u 
correspond to constant average energy for each mode. In 
Fig. [T] we plot the integral curves of the system defined 
in (fTO|) and (jlip for u = —vi and u = V2- 

For a feasible "bang-bang" strategy with only one in- 
termediate switching at t = ti, the appropriate control 
sequence is 



u{t) 



1, < = 

~vi, Q<t<ti 

V2, h < t < tl + t2 

l/7^ t = t'-p =h+t; 



(22) 



which is illustrated in Fig. 2(a) Applying the control 



boundary values in the opposite order does not transfer 
the state-space vector to the target. Note that the dis- 
continuities at the beginning and at the end of the pulse 
sequence are not implied by the maximum principle but 
from the initial and final conditions on the control u{t). 
If we ignore these boundary conditions and solve the cor- 
responding time-optimal problem, the minimum time ob- 
tained is a lower bound of the minimum time when these 
conditions are on. This bound is achieved with instanta- 
neous jumps of the control at the initial and final points. 
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0.5 1 1.5 2 2.5 
(a) Control function 



(b) Corresponding trajectory 



FIG. 2. The control function with one intermediate switching 
(panel a) and the corresponding trajectory (panel b) for vi = 
1,112 = 3 and 7 — 10. Dashed line corresponds to u = —Vi, 
solid line to u = 12. 



We next calculate the necessary time to reach the final 
point following the control strategy p2|) . Integrating 
and (dll) yields for t e [O.ti] 



xiit) = ./l + ^!l±lsinh\^,t), (23) 



while for t £ [ti,ti + t2 



x^{t) = J72 



sin^[y^(ti +i2 -<)]. (24) 



From (|2ip we find that the state-space equation of the 
first segment AB in Fig. |2(b)| is 



xi — WiX? H — 1 — V\, 

xi 



(25) 



since u = — fi and the starting point A(1,0) belongs to 
this segment. The corresponding equation for the second 



FIG. 3. The control functions with two intermediate switch- 
ings and the corresponding trajectories for iii = 1, W2 = 8, 7 = 
10. Dashed line corresponds to it = —v\, solid line to u = 12. 
Panels (a,b) show the intuitive solution with switching times 
given by (|32[) -(|34 [ ). while panels (c,d) show the optimal so- 
lution with at most two intermediate switchings where the 
switching times are calculated numerically. 



segment BF is 



xi 



1 V2 + 

T 



(26) 



since u = V2 and the final point F{j, 0) belongs to this 
segment. The two segments meet at point B{xf ,X2 )- 
Subtracting ^ from ^ we find that 



xf 



1 72^1 -I- 1 -f- 72(72^2 — 1) 



J^{vi+V2) 

Using the above value in and (|24l) we obtain 
1 



U = 



to 



: sinh 



1 



^"^{Vl + V2){vi -f 1) 

'w2(7^-l)(7^«i + l) 



/V2 V {Vl + V2){'j'^V2 — 1) 

The total necessary time is 



(1) 
/ 



tl+t2 



(27) 

(28) 
(29) 

(30) 



where the superscript denotes the number of intermediate 
switchings. 

We next show that when V2 is large enough we can find 
a control strategy with two intermediate switchings that 
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accomplishes the desired transfer in less time. Consider 
the following control sequence 



1, 

-Vl, 

1/7* 



t ^ 
0<t <ti 

ti <t <ti+t2 

ti+t2 <t <ti+t2+t3 



(31) 



where 



ti = - 



1 TT 



1 



■ sinh 



' {j^V2 - 1)(«2 +7'^Vl) 



(32) 
(33) 
(34) 



Time ii is chosen such that the first intermediate switch- 
ing takes place as close as possible to xi = (we explain 
later how this is related to minimizing the transfer time) , 
while t2 and are determined such that the target point 
F{j, 0) is reached at the final time. The control u{t) 
and the corresponding traje ctory for f i — l,V2 — 8 and 
7 = 10 are shown in Figs. 3(a) and |3(b)[ respectively. 
The total necessary time for this control policy is 



h+t2+ U 



Observe that for vi constant 



lim 



■ sinh 



^ 1^1(7^-1) 
vi + l 



lim 



(35) 

(36) 
(37) 



thus there is a value 



4(a)|we plot t^^^ 



V2 such that t^p 



In Fig 

Vl = \ and 7 = 



and t 



(2) 



<tV' for V2 > i'2- 



10. 



J aiiu uj as a function of V2 for 
Observe that for V2 > 6.786, the 



strategy including two intermediate switchings is faster. 

To understand intuitively why such a strategy can 
transfer the state vector to the final point in less time, 
we use the one-dimensional particle model where the po- 
sition xi and velocity X2 of a unit mass particle satisfy 
equations and pT|) . and refer to Fig. |3(b)[ If V2 is 
large enough, then the particle can be transferred rela- 
tively fast from starting point A, with position xi — 1 
and velocity X2 = 0, to point B, with < <C 1. At 
this point, the force term l/xf is very large and substan- 
tially accelerates the particle. When the particle passes 
through point C, with position xi — 1 same as the start- 
ing point A, it now has a nonzero velocity {x2 > 0) that 
allows it to travel faster at the final point F, with xi = 7. 
The repulsive potential at a;i = acts as a slingshot, re- 
sembling the gravitational slingshots used in aerospace 
engineering to alter the speed of a spacecraft. 

We emphasize that the values of ti , t2 and ^3 given by 
are not optimal but correspond to a suboptimal 



5p 
4.S\ 

■To 4 

3^ 

"3.5 
2.5 



5p 
4.5 
4 
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(a) Transfer times for one and (b) Optimal transfer time for 



two intermediate switcliings 



at most two intermediate 
switchings 



FIG. 4. Total transfer times for the presented strategies as a 
function of 112 for wi = 1 and 7 = 10. In panel (a) we plot 
of the strategy with one intermediate 

(2) 



the transfer time 



switching (solid line), as well as the transfer time of the 
intuitive solution with two intermediate switchings (dashed 
line). Observe that for 112 > 6.786 the second strategy is 
faster. In panel (b) we plot the transfer time t / for the optimal 
strategy with at most two intermediate switchings, retrieving 
similar results. 



intuitive solution. We can determine the optimal switch- 
ing times numerically if we vary ti (for specific ti, t2 and 
are automatically fixed from the requirement to reach 
the target point at the final time) and pick the value 
that minimizes the total transfer time. In Figs. 3(c) and 
|3(d)| we show the numerically calculated optimal solution 
with two intermediate switchings for the same parame- 
ter values that are used in Figs. 3(a) and |3(b)| for the 
intuitive solution. Observe that the two solutions are 
very similar but in the optimal one the first intermedi- 
ate jump takes pla ce slig htly earlier, before the Xi axis 
is reached. In Fig. |4(b)| we plot the total transfer time 
tf, obtained with this numerical method, as a function 



of V2 for Vl — 1 and 7 = 10. Comparing with Fig. 4(a) 
it is not hard to find that for 112 < 6.763 it is 

*(2) 



(2) 

while for V2 > 6.763 it is t/ w (actually tf < t 
as expected since the intuitive solution is suboptimal). 
In other words, for V2 < 6.763 the numerically calculated 
optimal solution has only one intermediate switching (the 
method gives ti = indeed), while for V2 > 6.763 is very 
close to the intuitive solution given by (l5T]) - p4l) . Note 
(2) 

that since tf < tj, the transition value of V2 just found 
(6.763) is slightly lower than the one found above (6.786) 
using the intuitive solution. 

For larger values of V2 it is possible to find strategies 
with more than two intermediate switchings with shorter 
transfer times, depending on the value of the target co- 
ordinate 7. For example, consider the strategy with 2n 
intermediate switchings whose corresponding trajectory 
is shown in Fig 5(a) It is actually composed by n seg- 
ments with two switchings, see Fig. |5(b)[ The necessary 
time to travel the segment starting from (/3i_i,0) and 
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(a) Trajectory with 2n 
switchings 



(b) Segment with two 
switchings 




(c) Transfer time t'^"' as 



V2 



f 

oo for optimal /3i 



40 60 

(d) Exact transfer times for 
two, four and six switchings 



FIG. 5. A trajectory with 2n intermediate switchings (panel 
a) composed by n segments with two switchings (panel b). 
Dashed line corresponds to u = —vi, solid line to it — V2- 
Note that /3o = 1 and /?„ = 7. In the limit ?;2 — ^ 00 the 
minimum transfer time corresponding to the optimal choice 
of Pi is plotted (panel c) for 7 = 10 (solid line) and 7 = 50 
(dashed line) as a function of n. Observe that for larger values 
of 7 the minimum is shifted towards larger values of n. For 
7 = 10 it is achieved for n = 2, i.e. four intermediate switch- 
ings. The exact transfer times t^p , t j'-' and , as calculated 
from (|38|) (valid at all scales of V2) are plotted as functions of 
V2 for 7 = 10 and vi — 1. Observe that for 112 > 43.32 the 
four switchings strategy (dashed line) becomes indeed optimal 
among the control policies that we consider. 



ending at {(3i, 0) is 



if (ft_i,A) = ^1 



1 



Pi-Iy/V2 



1 



2v^ 



(38) 



where 

i2(a,/3) 



1 lv,iP^-a^)ia^(3^V2~l) 



(3^ivi+V2){a'^vi + l) ' 



In the limit U2 we obtain 




The total transfer time for the strategy with 2n switch- 
ings is 



(40) 



where /3o = 1 and /3„ — 7. We can find the opti- 
mal /3i,i = l,2,...n — 1 that minimize t^^"^ using dy- 
namic programming. Suppose that we know the optimal 

/3i, i = 1, 2, ... n — 2 and we want to find Pn-i- This 

(2) 

variable appears only in the terms (/3„_2, /3n-i) and 

tf\l3n-i,l3n) of the sum (gO]). Using Eq. ^ to approx- 
imate these terms and equating with zero the derivative 
of their sum with respect to this variable, we find that 
the optimal /3„_i satisfies P^-i = Pn-2Pn in the limit 
f2 — 00. It corresponds to a minimum since the second 
derivative can be easily found to be positive. Working 
analogously we get 



for i — 1, 2, . . . n - 

/3, =7'/",z = l,2, 



transfer time t 



,{2n) 



(2n) 
/ 



/32 = (41) 

l,/3„ = 7, we obtain 
The minimum value t i^"^ of the 



1. Since /3o 
. . n 

as V2 — > 00, is 



- + ^72/" - 1 + sin 



, (42) 



where note that each of the n segments is traveled in 



equal time. In Fig |5(c)| we plot t^j"^ in units of 
as a function of n for 7 = 10 (solid line) and 7 = 50 
(dashed line). Observe that for larger 7 the minimum of 

^fm shifted towards larger values of n, i.e. the particle 
passes more times close to a:i = to acquire more speed 
and thus reach faster the more distant target point. For 
7 = 10 the minimum is obtained for n ~ 2, i.e. a four 
switchings strategy. In Fig |5(d)| we plot t^j"\n = 1, 2, 3, 
for 7 = 10 and vi = 1, using the exact formula ([55]) and 
not the approximation Observe that for V2 > 43.32 
the four switchings strategy is indeed the optimal among 
the control policies that we considered. 

Although the time-optimal control has "bang-bang" 
form, as shown above, such discontinuous changes in u{t) 
are unrealistic and difficult to implement experimentally. 
To overcome this problem, in the next section we use a 
powerful numerical optimization method that allows us 
to find realistic time-optimal solutions, following the path 
introduced in 



IV. REALISTIC SOLUTIONS USING A 
PSEUDOSPECTRAL NUMERICAL METHOD 

Pseudospectral methods were developed to solve par- 
tial differential equations and recently adapted to solve 
optimal control problems, see for example |27H3ll | and our 
recent work for optimal pulse design in Nuclear Magnetic 
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-0.5 0.5 

(a) Uniform Grid 




-0.5 0.5 

(b) LGL Grid 



FIG. 6. The A*' = 16 order interpolation of the function 
f{x) = l/(16x^ + 1) based on a uniform grid (panel a) 
demonstrates the Runge phenomenon whereas the interpo- 
lation based on the LGL grid (panel b) does not. 



Resonance spectroscopy [20]. They are used to convert 
a continuous-time optimal control problem to a discrete 
nonlinear programming problem, which can be solved by 
many well developed computational algorithms. 

These methods involve the approximation of the con- 
trol and state functions, u(t) and x(t), by orthogonal 
polynomial basis functions on the domain [—1, 1]. Using 
such a basis leads to spectral accuracy, namely, the /c*^ 
coefficient of the expansion decays faster than any inverse 
power of k [35| , permitting the use of relatively low order 
polynomials in the approximations. 

In order to apply such a method, the first step is to 
transform the optimal control problem from the time do- 
main t e [0;tf] to r e [—1,1] using the simple afffne 
transformation r(t) = (2t ~ tf)/tf. In a redundant use 
of notation, we make this transition and reuse the same 
time variable t. The system equations become 

• k 

2 



■X2 



-uxi H o 



(43) 
(44) 



with boundary conditions (a;i(— 1), a;2(— 1), u(— 1)) — 
(1,0,1) and (xi(l),X2(l),w(l)) = (7,0,1/7"). 

According to the Chebyshev Equioscillation Theorem 
[33I ] the best A^*^ order approximating polynomial to a 
continuous function on the interval [—1, 1], as evaluated 
by the uniform norm, is an interpolating polynomial. 
Since any N*"^ order interpolating polynomial can be rep- 
resented in terms of the Lagrange polynomials, we use 
these functions to express the interpolating approxima- 
tions of the state and control functions, x(t) and u{t). 
Given a grid of -I- 1 interpolation nodes within [—1,1], 
r = {to < ti < ■ ■ ■ < tN}, the Lagrange polynomials 
{£k}, k G {0, 1, . . . , A^}, are constructed by 



N 



- n 



i = 



{t - t,) 
{tk — ti) 



They form an orthogonal basis with respect to the 
discrete inner product {p,q) — '^^-Qp{tk)Q{tk), while 
^k{ti) = Ski holds at the grid nodes [34i] . 




1 2 3 4 5 6 



(a) Optimal control, M = 00 (b) Corresponding trajectory 




(c) Realistic control, M 



10 (d) Corresponding trajectory 



FIG. 7. Control functions calculated by the pseudospectral 
method for the same parameters as in Fig. [2] without (M = 
00, panel a) and with (M — 10, panel c) slope restriction. 
The latter case requires a larger transfer time, as expected. 
The corresponding trajectories are also shown (panels b,d). 



For an arbitrary selection of nodes, as the order of ap- 
proximation N gets large, Runge oscillations near the 
endpoints of the [—1, 1] domain may occur |35;j, as shown 
in Fig. [S] In order to suppress this phenomenon and 
increase the accuracy of the approximation, we use the 
Legendre-Gauss-Lobatto (LGL) nodes, which are the end 
points ^0 = — 1, ^Af — 1 and the roots of L]^{t), the deriva- 
tive of the A^**^ order Legendre polynomial 36] . The cor- 
responding grid is r'^'^'^ = {ti : to = -l,iAr(t)]t^ = 
0, i = 1, . . . A^ — 1, ^AT = 1}. In this case the Lagrange 
polynomials £kit) can be expressed as 



4(i) 



{t^-l)LNit) 



N{N + l)LN{tk) t-tk 



(45) 



where {tk} e T^'^^ , fc = 0, 1, . . . , A^. 

The N^^^ order interpolating approximations of the 
state trajectory and control functions with respect to the 
same grid are. 



At) 

u{t) 



iNy^it) = J2k=0^kikit), 



lNu{t) = Y.k=0^kik{t), 



(46) 

(47) 



where and Uk are not only the coefScients of the ex- 
pansions, but also the function values at the fc*^ node 
due to the definition of the Lagrange polynomials [2?} . 
From the interpolation as in (|46p. we have 



N 



dt 



/jvx(t) -^Xfe4(t). 



fe=0 



8 



10 




(a) Optimal control, M = oo (b) Corresponding trajectory 



10 




(c) Realistic control, M = 10 (d) Corresponding trajectory 



FIG. 8. Control functions calculated by the pseudospectral 
method for the same parameters as in Fig. [3] without (A/ — 
oo, panel a) and with (Af — 10, panel c) slope restriction. 
The latter case requires a larger transfer time, as expected. 
The corresponding trajectories are also shown (panels b,d). 



lowing dynamic constraints 

DikXik = y X2i (50) 

fc=0 

^ DikX2k = ^ i-UiXu + -\-] (51) 

fc=o . ^ ^i,*^ 

for i = 0,1, . . . , N , with Xj = {xii,X2i) ■ To prevent 

unrealistic discontinuities in u{t), we impose the following 

slope restriction 

— < M (52) 

U+i — ti 

for i = 0,1, . . . , N — 1, where M characterizes the max- 
imum allowed slope of the control function. The cor- 
responding finite-dimensional constrained minimization 
problem is to find minimum tf and {ui} with —vi < Ui < 
V2, such that the above algebraic relations and the bound- 
ary conditions {xio,X2o,ua) = (1,0,1), ixiN,X2N,UN) = 
(7, 0, 1/7^) are satisfied. Solvers for this type of problems 
are readily available. In Figs. [7]and|8]we plot the optimal 
controls and the corresponding trajectories calculated by 
the pseudospectral method, for the same parameter val- 
ues as in Figs. [2] and |3l respectively, with and without 
slope restriction. 

V. CONCLUSION 



Using (j45|) and special recursive identities for the deriva- 
tive of Legendre polynomials [3l|, we have at the LGL 
nodes U £ T^gl^ i = 0,l,...,N, 



N 



N 



dt 



(48) 



fc=0 



*:=0 



where Dik are ifc"^ elements of the constant (N + 1) 
{N + 1) differentiation matrix D defined by [SJl 



/ -Ljv(t.) 



A, 



LN(tk) ti-tk * ^ ^ 

i^k^O 
i = k = N 
otherwise. 



N{N+1) 
4 

N{N+1) 
4 



(49) 



The pseudospectral method is a collocation method in 
which the state dynamics is enforced at the LGL nodes. 
Using dig), dm), (gni), (gZl) and (gS]), we obtain the fol- 



In this paper we used optimal control theory to show 
that minimum time frictionless atom cooling in harmonic 
traps is achieved when the trap frequency changes in a 
"bang-bang" manner, even if the trap is allowed to be- 
come transiently an expulsive parabolic potential. Us- 
ing this fact we calculated estimates of minimum cooling 
times for control strategies with various numbers of fre- 
quency jumps. Finally, we employed a pseudospectral 
optimization method to find realistic solutions without 
discontinuities, appropriate for experimental implemen- 
tation. The above results and techniques are not re- 
stricted to atom cooling but are applicable_ to areas as 
diverse as adiabatic quantum computing [3| and finite 
time thermodynamic processes [3]. 
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